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ABSTRACT 



This work examines the three-dimensional dynamic response 
of a nonlinear fast reactor with temperature-dependent feed- 
back and delayed neutrons when subjected to uniform and local 
disturbances. The finite element method was employed to re- 
duce the partial differential reactor equation to a system of 
ordinary differential equations which can be numerically in- 
tegrated. A program for the numerical solution of large 
sparse systems of stiff differential equations developed by 
Franke and based on Gear's method solved the reduced neutron 
dynamics equation. Although a study of convergence by refin- 
ing element mesh sizes was not carried out, the crude fi- 
nite element mesh utilized yielded the correct trend of neu- 
tron dynamic behavior. 
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I. INTRODUCTION 



The nuclear reactor with delayed neutrons and tempera- 
tu re-dependent feedback is a nonlinear system, whose response 
to uniform and local disturbances differs greatly from that 
of a linear reactor. The prompt temperature feedback model 
and one average group of delayed neutrons were incorporated 
in this work. Practically all neutron dynamics analysis to- 
day deals with two-dimensional geometry, implying a symmetry 
in neutron dynamics behavior during transient [1,2]. This 
symmetry assumption, however, could be unrealistic in the 
safety analysis of nuclear reactors. A unique feature of 
this work is the consideration of the three-dimensional depen- 
dence of the neutron flux. No symmetry assumptions were im- 
posed on the problem. This work investigated neutron dynamics 
under uniform and local disturbances in the core. A uniform 
initial flux throughout the interior of the reactor was im- 
posed. The finite element method (FEM) was employed to solve 
this nonlinear i ni ti al -boundary value problem. The FEM is 
quite effective in handling discontinuous forcing functions 
thereby making it particularly suited for examining the ef- 
fects of localized pertu rba ti ons and space-dependent feedbacks. 
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1 1 . THE NUCLEAR REACTOR WITH TEMPERATURE 

DEPENDENT FEEDBACK AND DELAYED NEUTRONS 



A. THE PHYSICAL SYSTEM 

The system under consideration is a fast reactor of cy- 
lindrical geometry that is composed of two different regions. 
The reactor core or region I is cylinderical in shape and is 
fueled by U-235. Region II or the reflector region complete- 
ly surrounds the core and is composed of U-238. Both regions 
were assumed to be homogeneous. Table I lists the physical 
properties and geometry for each region. A schematic of the 
fast reactor geomenry is shown in figure 1. 

Temperature and delayed neutron effects were taken into 
account. Also, a one-velocity or one-group model was assumed, 
thereby making the velocity independent of spatial or tem- 
poral effects. The delayed neutrons were considered only in 
region I since region II was assumed to be a non-multiplying 
medium. The temperature effects were also assumed to be only 
in the core region. 

In general form, the one-velocity neutron diffusion equa- 
tion is [ 3 ] 

= vDV^N - Z vN + S(r,t) (1) 

o L a 

where r = (x,y,z) 

D = neutron diffusion coefficient 
Z = macroscopic absorption cross-section 

d 

N(f,t)dV = number of neutrons in a volume element dV at 
a point r at time t 
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TABLE I 



Physical Constants 



SYMBOL 


DEFINITION 


VALUE 




radius of Region I 


60 cm 




total reactor rajdius 


90 cm 


% 


height of Region I 


160 cm 


hi 


total reactor height 


220 cm 


V 


neutron speed 


7 

4 . 8x10 cm/ sec 




core neutron diffusion coefficient 


0.913 cm 


“ll 


reflector neutron diffusion 
coefficient 


1,200 cm 




core neutron absorption cross 
section 


0.01401 cm"^ 


1 

“II 


reflector neutron absorption 
cross section 


0.0040 cm ^ 


V 


number of neutrons per fission 


2.54 


•U 

h 


critical fission cross section 


0.005736 cm"^ 


6 


delayed neutron fraction; 
dollar reactivity 


0.00642 


£ 


fission energy 


7.652x10"^^ 

cal/fission 


h(^) 


modified convection heat 
transfer coefficient 


0.0632 . 
cal/ (cm^ sec °C) 


a 


tenperature coefficient 


-0.004/°C 


1 


abundance-weighted mean decay 
constant 


0.4350 sec ^ 



I 

I I 



Y 




core 
ref 1 ector 




Figure 1. Schematic of cylindrical reactor. 
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vDV^NdV = number of neutrons diffusing into dV per 
unit time at time t 

E^vN dV = number of neutrons absorbed in dV per 
unit time at time t 

S(f,t)dV = number of neutrons produced in dV per 
unit time at time t 



The neutron number density is related to the neutron flux by 
the expression [4] 



(J)(r,t) = vN(r,t) (2) 

where (|)(r,t) is the flux at time t . The neutron diffu- 
sion equation in terms of the flux is depicted by 

i = DV"<1) - Eg(j) + S(F,t) (3) 

B. PROMPT AND DELAYED NEUTRONS 

The source or production term in equation (3) is composed 
of the contributions of the prompt and delayed neutrons. The 
majority of the fission neutrons are prompt neutrons that 
appear almost instantaneously (within 10"^ second) on fission. 
Assuming a fast neutron non-leakage probability of unity, the 
prompt neutron source is described by 

Sp(r,t) = (l-B)K^(r,t) <l)(r,t) (4) 

where 

“ infinite medium multiplication factor 
S = total fraction of delayed neutrons 
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The fraction of delayed neutrons is very small (note 

that 6 = 0.00642 from Table I). However, they have a very 

significant effect on the reactivity because their mean life 

times are long. Without delayed neutrons, reactor control 

would not be possible. These delayed neutrons are born in 

the decay by neutron emission of nuclei produced following 

the 6-decay of certain fission fragments. For example, the 

87 86 

6-decay of the fission fragment Br leads to Kr plus a 

8 7 

neutron. Nuclei such as Br whose production in fission 
eventually leads to the emission of a delayed neutron are 
called delayed neutron precursors [4]. 

There are six main groups of delayed neutrons. Each 
group is classified according to its decay constant. The de 
layed neutron source term is portrayed by 

6 

Sj(r,t) C^. (f,t) (5) 

i = l 



where 

t h 

= decay constant of the i group 
C^.(r,t) = density of the i^^ precursor 

Assuming that the fission fragments do not migrate appre 
ciable distances and assuming a non-circulating fuel reactor 
[3], the precursor density is delineated by 



9C^. (r ,t) 
Tt~ 



6.K Z 4> - X.C. 
1 " a 11 



( 6 ) 



where 6^- 



i s 



the fraction of delayed neutrons of the i 



th 
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group. The solution to the precursor equation is in terms 

of a time integral expressed by 

t 

C.j(r,t) = y* e"^i^^"^ ^ . t )(}) ( r , t ) dt ' (7) 

0 

Inserting equation (7) in equation (5) yields the delayed 
neutron production term as 



=i; 

1=1 






/ 



g-^-j (t-f ) 



K„(r,t)(})(f,t)dt' 



° (8) 

For convenience the six main groups of delayed neutrons 
were considered as one group. This was accomplished by us- 
ing the abundance-weighted mean decay constant defined by 



A 




i = l 



6 .A . 
1 1 



(9) 



This is a reasonable approximation since many of the impor- 
tant phenomena of nuclear reactor dynamics can be character- 
ized satisfactorily by combining all the emitters or precur- 
sors into one, two, or, at most, three effective groups [3]. 
Replacing A^. with the abundance-weighted mean decay con- 
stant showed the delayed neutron source term to be 

t 

S^(?,t) = 6AZ^ y' ^K^(r,t)(J)(r,t)dt' (10) 

0 

C. DOPPLER EFFECT AND TEMPERATURE FEEDBACK MODEL 

The Doppler effect in fast reactors is due to the tempera- 
ture broadening of many closely spaced high-energy resonances 
in both the fission and parasi tic-absorption cross-sections. 
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These resonances basically mean that as the temperature in- 
creases, the number of neutrons that are absorbed increases 
and, similarly, the number of fission events increases. 
These nonproductive and productive processes compete in a 
complicated manner, and the net effect may be either an in- 
crease or decrease in reactivity [3]. For most fast reac- 
tors, the effect is negative, and it will be so assumed in 
this work. 

The reactivity change with respect to the fuel tempera- 
ture change is modeled by [2] 



where a, b, and c are parameters determined from experimen- 
tal or neutronic calculations, m is an integer, and T is 
the current fuel temperature. For most fast reactor cores 
(ceramic-fueled cores), TdK^/dT is very close to being con- 
stant over a wide range of temperatures. Therefore, it is 
assumed that a and c are zero and the Doppler coefficient 
is expressed by 



dT 



CO 



aT'^/^ + bT“^ + cT^*""^ ^ 



( 11 ) 



dK 



oo 



( 12 ) 



The following initial conditions were used: 




(12a) 

(12b) 
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The reactivity model is then expressed by 



K.(f.t) - K” + b Hn (J-) (13) 

0 



where 

K° = multiplication factor at steady-state 
V = number of neutrons produced per fission 
Tq = original fuel temperatures 
b = Doppler constant 

The definition of K® is 



K 



O 

cx> 




(14) 



where 

★ 

= critical fission cross-section 

c 

Z = macroscopic absorption cross-section of the 
core 



The rise in fuel temperature is depicted by 

e(r,t) = T(r,t) - T^(r) (15) 

This is also described by the integral [5] 

t 

9(r,t) = y* f(t-t') ij;(r,t)dt' (16) 

0 

where f(t-t') is the feedback kernel and is dependent upon 
the type of temperature feedback model used. i()(r,t) is the 
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rise in neutron flux above the steady state and is expressed 
by 



ii^(r,t) = 4>(f,t) - 



(17) 



where 1 s the neutron flux at steady state. 

There are three types of temperature feedback models that 
can be considered here. The first is known as "Newton's 
feedback" which determines the reactor temperature by Newton's 
law of cooling. The second temperature feedback model is 
called the adiabatic feedback model. This represents the 
temperature for the loss of coolant case. The third model is 
the prompt temperature feedback model in which the fuel tem- 
perature follows the behavior of the neutron flux without de- 
lay [5,6]. The prompt feedback model was employed in this 
analysis. The feedback kernel for this temperature feedback 
model is 



where K is an energy production operator with units of °C 
per unit flux and y » a dimensionless quantity, is related 
to the mean time for heat transfer to coolant. 

Inserting equation (18) in equation (16) and performing 
the integration produced a rise in temperature of 



f(t-t' ) = ^ 6(t-f ) 



(18) 



9(r,t) = tp(r,t) 



(19) 
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From equation (19) the temperature of the fuel is 

T(r,t) = ^ i|^(r,t) + T^(?) (20) 

The ratio of the current temperature to the steady -state tem- 
perature is therefore 

r- = t + 1 (21) 

In this work was considered to be constant throughout 

the core. Incorporating equation (21) into equation (13) 
gave the reactivity model of 

K^(r,t) = K° + b Jin [;j^ ,t) + 1] (22) 



D. FIELD EQUATIONS 



Before establishing the field equations it is desirable 
to express equation (3) in terms of the rise in flux. Using 
equation (17) in equation (3) and grouping terms yielded the 
following diffusion equation; 



1 ^ = 
V 3t 



+ 



+(1-6)K ^ 

a ^ a a 









0 



_ 



(23) 



The second bracketed term in equation (23) is identically 
equal to zero since it is the steady state portion of the 
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diffusion equation. The rise in neutron flux above its 
steady state value is therefore expressed, for the core, by 

7 It " 

+ eXz^ ^ 24 ) 

0 

and, for the reflector region, by 

H! ' - h* <25) 

Inserting the reactivity model into equation (24) yielded 



i|l = * (1-B) 

t 

+ (l-B)bZ^ {jln[^ + 1]} ip + SAK° y* e'^ ^ ^ ' ^ij;dt ' 



+ BXZgb { y* 
0 



)[Jln(^ + 1)] ipdt'} 

^ 0 



(26) 



For the reflector, the last four terms of equation (26) are 
zero. The effects of the temperature on the delayed neutrons 
were neglected in this work. The field equations can now be 
expressed, for the core, as 



- vDV^ij; + [vZ^ - v(l-3) vZ^]ijj 
t [-(I-B)vs b] [ln(^ + 1)]* 

+ [-BXvZ^v] [ y'%~^(^’^'^ij;df] = 0 (27) 
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and, for the reflector, as 



II - vDV^ij; + = 0 (28) 

The non-linear terras of the core field equation will be 
linearized accordingly. In raore compact form, equations (27) 
and (28) became 

II - civ^ij; + c2^|, + c4[Jln(^ +1 )]ij^+c5[/^ e’^ ^ ' ^4^dt ' ] = 0 

(29) 

and 

II - clV^4; + c2i|; = 0 (30) 

where the meanings of the coefficients cl, c2, etc., are 
obvious for the core and reflector. 

Equations (29) and (30) were subjected to the following 
condi ti ons : 

boundary condition: 

4'R(rg,t) = 0 (31) 

where are coordinates of points on the outer surface of 

the reflector and the subscript R refers to the reflector 
continuity of flux: 

^g(r,t) = ,t) (32) 

where r^ are coordinates of points on the core- ref 1 ector 
interface and the subscript c refers to the core. 
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III. FINITE ELEMENT 



A. INTRODUCTION 

The application of the finite elements of nonlinear con- 
tinua has been mostly in the field of solid mechanics. 

Prior work [1] using the finite elements of nonlinear con- 
tinua on a nonlinear reactor dynamics problem has been suc- 
cessful. The FEM in this work was utilized to reduce a non- 
linear partial differential equation of the nuclear reactor 
to a system of nonlinear ordinary differential equations in 
time. The time integration was accomplished by using a com- 
puter program for the numerical solution of stiff differen- 
tial equations developed by Franke [7]. In order to minimize 
computer storage requirements, an optimum compacting scheme 
(DCS), described by Ref. 8, was adopted. 

The finite element models of operator equations are gen- 
erally classified into three categories: a) the variational 
finite element models such as the Ritz method, b) the weigh- 
ted residuals method such as the method of Galerkin, and c) 
the direct finite element models which are not based on 
functional minimization. From experience in structural 
mechanics, the most effective method for generating accept- 
able finite element models of nonlinear equations is the 
Galerkin method [1]. This work adopted the method of Galer- 
kin in a finite element approximation over the spatial 
domain of the field equations. 
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B. THE METHOD OF GALERKIN 

The Galerkin method is a special case of the method of 
weighted residuals. It involves a rational choice of 
weighting function that is consistent with the type of fi- 
nite element approximation considered. Indeed, the weight- 
ing functions chosen are the basis or shape functions em- 
ployed in the finite element approximation. There are two 
favorable characteristics of the Galerkin method which 
makes it attractive. The first attribute is its amenability 
to integration by parts. This supplied the freedom of using 
a lower order finite element than might be otherwise possibl 
The second favorable characteristic of the Galerkin method 
is that the symmetric operators in the field equations trans 
form into symmetric matrix operators. Both these attributes 
are attractive for computational purposes. 

Consider the ini ti al -boundary-val ue problem 

^ = 1 ^ - f(r,t) (33) 

where contains the nonlinear operators. According to 

the spatial finite element discretization, the solution of 
equation (33) is in the form of an union of an N-term approx 
imation given by 

R 

>l^(r,t) = i(r,t) = LJ iljJt) G,(f) , J = 1,2,...,N (34) 

J = 1 

where N is the number of system degrees of freedom (i.e., 
number of coordinates), and Gj(r) are the system or global 
basis functions which span the space of the approximate 
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solution ii)(r,t) [1]. The Einstein summation is used. The 
global basic functions are "pyramid functions", each of 
which has a prescribed functional description over a sub- 
domain of the system and is zero elsewhere [9]. The unknown 
coordinate functions t|jj(t) are the time-dependent magni- 
tudes of the approximated flux ijj(r,t) and/or its deriva- 
tives at discrete nodal points [10]. 

The residual function, R(r,t) , is defined such that it 
is identical to zero when ijj(r,t) is equal to the exact solu 
tion. The residual function is expressed by 



The Galerkin orthogonality condition (using the basis func- 
tions as weight functions), when applied to the residual func 
tion, requires that 



Vol 

From the f i el d equati ons , the residual function for the core 
i s 



R(r.t) = II -iij; - f 



(35) 




(36) 



R(f,t)= II - clV^ij; + c2ip + c4[2,n(|| + l)]>j; 




(37) 



0 



and for the reflector 



R(r,t) = Iy - clv^ij) + c2i 



(38) 
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Using equation (34) and applying Galerkin's orthogonality 
condition produced a core equation of 

/6i{GjiIjj(t)-V2Gj(cl-(|;)j + Gj(c2-ii;)j 

Vol 

A . 

+ [ / hc5-i|))j dt']Gj} = 0 (39) 

0 

where I = 1 , 2 , . . . , N 
J = 1,2,. ..,N 
K = 1 ,2, . . . ,N 

The reflector has a similar equation without the nonlinearity 
and is expressed by 

j Gj{Gji^-j(t) - v2Gj(cl-(^)j + Gj(c2-tp)j} dVol = 0 (40) 

Vol 

C. THE ELEMENT 

A three-dimensional quadratic isoparametric element was 
employed in this work. The parent element is a triangular 
prism or solid wedge with straight sides. The element shape 
functions are expressed in terms of area coordinates in the 
plane of the triangle and by an isoparametric coordinate 
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along the prism axis. Figure 2 shows the parent element. 
This element was chosen because of the ease with which it 
fits the cylindrical structure when it is transformed into 
a curved element. This type of element has been used before 
as filler elements [11]. 

The area coordinates are defined by area ratios. Con- 
sider the triangle shown in figure 3. An arbitrary point P 
within the triangle defines three subareas designated by , 

and A^. The ratio of each of the subareas to the total 
area is known as an area coordinate. In equation form, the 
area coordinates are 



where A is total area of the triangle. L.| , L 2 , and 
are the natural coordinates for a triangle. The requirement 
that the sum of the subareas be equal to the total area is 
obviously satisfied by the identity 



In the plane of the triangle, only two of the area coordi- 
nates are independent. 

The isoparametric coordinates are best visualized by con- 
sidering the rectangular prism shown in figure 4. Isopara- 
metric coordinates are normalized coordinates such that 
their values on the faces of the rectangle are + 1. The Cn? 
axes are in general not orthogonal. They are orthogonal 



L2 = k^/k 
L3 = A3/A 



L, = A, /A 



(41a) 

(41b) 

(41c) 




(42) 
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Figure 2. Quadratic triangular prism parent element 
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Y 




Figure 3, Definition of area coordinates 
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? = 1 




Figure 4. Isoparametric coordinates 
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only in the special case of a rectangular prism element 
[12]. The element basis functions use the c coordinate 



in the prism axis and the and coordinates in 

the plane of the triangle. 

The parent element, as shown in figure 2, has 15 local 
nodal points around its periphery. As such, there are 15 
basis functions which are given below [11]: 

Corner nodes: (nodes 1, 3, 5, 10, 12, 14) 

= 1 L^(2L^-l)(l+0 - i L^(l-cM 
N3 = \ L2(2L2-1)(1+c) - j L2(1-cM 
= 1 L 3 ( 2 L 3 - 1 )( 1 +c) - j L 3 O-CM 
N^q= \ L^(2L^-l)(l-c) - ^ L^(l.-<;2) 

N^2= i L2(2L2-1)(1-?) - ^ 12 ( 1 - 4 ^ 

1 L3(2L3-1)(1-c) - I L3(1-cM 

Midside nodes of rectangles: (nodes 7, 8, 9) 

N 7 = L^(l-cM 
Ng = 12(1-4") 

Ng = L 3 ( 1 - 4 ") 

Midside nodes of triangles: (nodes 2, 4, 6, 11, 13, 15) 

N 2 = 21 ^ 12 ( 1 + 4 ) 

N4 = 21213(1+4) 

Ng = 2131^(1+4) 

21 ^ 12 ( 1 - 4 ) 

*^13^ ^^2''3^ ^ 

Ni 5 = 2131 ^( 1 - 4 ) 



31 



The coordinates of each local- node, in terms of L^, 

and €, are listed in table II. These element basis functions 

< N > define the geometry of the element. Note that they 
satisfy the relationship 

1 , at node i 

0 , at other nodes (43) 

On the element level, the variation of the unknown func- 
tion \p is approximated by 

= < N •> {ij; }® (44) 

where ■ < N'> = row vector of element shape functions 

0 

{ } = column vector of time-dependent nodal 

values of ip® 

To satisfy continuity requ i rements , the shape functions 

< N‘> have to be such that the continuity of the unknown 
function ip is preserved in the parent coordinates [11]. 

The shape functions < N > which characterize the element 
geometry and the shape functions < N'> which describe the 
unknown function do not necessarily have to be the same. 
There is no requirement that the nodal values be associated 
with the same nodes which were used to define the element 
geometry, though in practice it is often the case. Consider 
for example, the illustrations in figure 5, If the nodes 
defining the element geometry and the nodes defining the 
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TABLE II 



Coordinates of Local Nodal Points 



Local Node 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 



1 

0 

0 

0 

% 

1 

0 

0 

1 

h 

0 

0 

0 

% 



0 

1 

0 

0 

0 

1 

0 

0 

% 

1 

0 

0 



5 



0 

0 

0 

1 

0 

0 

1 

0 

0 

0 

1 

3s 



1 

1 

1 

1 

1 

1 

0 

0 

0 

-1 

-1 

-1 

-1 

-1 

-1 
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(a) Isoparametric 




• - geometry nodes 
O - variable function nodes 

Figure 5. Element classification 



3A 



unknown function are identical, the element is known as an 
isoparametric element. This means that the shape functions 
describing the geometry and the shape functions describing 
the unknown function Tp are equal, or 

< N’> = < N > (45) 

If there are more nodes defining the geometry than nodes de- 
fining the variable function, the element is called a super- 
parametric element. Using more nodal points to define the 
unknown function than to describe the geometry of the element 
results in a subparametric element [11]. This work utilized 
the isoparametric element classification. 

D. DIVISION OF THE SYSTEM INTO ELEMENTS 

In three-dimensional space, the division of the sys- 
tem into discrete elements is difficult to visualize. It is 
virtually impossible to show every nodal point of the system 
in one schematic. To present a clear view of the discretized 
domain, a "layer" approach was adopted. The first finite 
element grid or mesh employed here consisted of 128 elements. 
Under this grid (mesh I), the reactor was divided into four 
layers as shown in figure 6. Each layer was composed of 32 
elements. The first and fourth layers were each 30 cm in 
height and each contained entirely reflector elements. The 
second and third layers were each 80 cm in height and together 
they encompassed the entire core plus the remaining reflector 
elements. Each layer, in turn, was partitioned into three 
horizontal (xy) planes. The top plane included all the global 
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z 




I - core II - reflector 

igure 6. Layers of mesh I 



110 



80 



=0 



=-80 



-110 
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or system nodes corresponding to local or element nodes 1 
through 6. The middle plane contained all the system nodes 
corresponding to the local nodes 7, 8, and 9. System nodes 
corresponding to element nodal points 10 through 15 comprised 
the bottom plane. To fix ideas, the three nodal planes of 
the first layer of mesh I are shown in figures 7, 8, and 9. 

In this work there was only one curved side which was in 
the plane of the triangle, as shown in figure 10. The first, 
seventh, and tenth element nodes were each arbitrarily assigned 
to have as an opposite side the curved side of the triangle. 

The remaining local nodes in each of the respective planes 
of the element were then numbered consecutively in the counter- 
clockwise direction. 

At this point it is appropriate to define connectivity. 

The connectivity of an element is a row vector array that 
relates the local nodes to system nodal points. The connec- 
tivity lists the system nodal points that are "connected" to 
form the element domain in the sequence of local nodal number- 
ing. Thus, the connectivity matrix for mesh I is a matrix 
of size 128 x 15. To illustrate, the connectivity of element 
number 106 is 

<235, 210, 193, 194, 195, 211 , 266, 176, 177, 283, 152, 10, 11 , 12, 53> 

A second finite element mesh (mesh II) consisting of 192 
elements was developed. It contained the same number of ele- 
ments per layer as mesh I. However, mesh II has six layers. 

The second and third layers of mesh I were each divided in 
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Number - element number 



Figure 7. Top nodal plane of the first layer of mesh I, 
z = 110 cm 
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265 




Number - element number 



Figure 8. Middle nodal plane of the first layer of 
mesh I , z = 95 cm 
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Number - element number 



Figure 9. 



Bottom nodal plane of the first layer of 
mesh I , z = 80 cm 
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5 





Figure 10. Local nodal numbering of curved 
el ements 
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half, thus forming the two additional layers of mesh II. 

The connectivity matrix and the coordinates of each global 
node for mesh I are given in Appendix A. Appendix B lists 
the connectivity matrix and nodal coordinates of mesh II. 

The reader can construct the different layers and elements 
without much difficulty by using the nodal coordinates and 
the connectivity matrix. A mesh generator was not utilized 
in this work. 

E. COORDINATE TRANSFORMATION 

The use of a quadratic or higher order element permits 
the transformation or mapping of the straight-sided parent 
element into an element with curved sides. Distorted or 
curved elements provide a better fit to curved domains than 
linear elements, and thus a smaller number of elements is 
required to represent the structure adequately. 

The transformation from cartesian coordinates to curvi- 
linear coordinates can be accomplished by employing a one-to- 
one correspondence defined by [11] 



X 

y 


= f 


n 


or f ( 




(46) 


z 




c 

1 - 









The element shape functions are utilized to achieve this 
transformation via the relation 
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0 

where n is the number of element nodes, and the element 
shape functions N^. are in terms of local coordinates. Each 
set of local coordinates corresponds to only one set of car- 
tesian coordinates. 

In performing the transformation, the compatibility re- 
quirement must be met. The transformation into the new, 
curved elements should leave no gaps between adjacent elements. 
If two adjacent elements are generated from parents in which 
the element shape functions satisfy continuity requirements, 
then the curved elements will be contiguous [11]. For the 
isoparametric element, uniqueness of coordinates ensures com- 
patibility. Continuity is assured when adjacent elements 
are given the same sets of coordinates at common nodes. 

Since the element shape functions are in terms of local 
coordinates, an element of volume, dxdydz, must be transformed 
into an element of volume expressed in local coordinates. 

This is achieved through the use of the Jacobian matrix de- 
fined below. Using the chain rule, the relationship between 
^,n,? and a corresponding set of cartesian coordinates 
x,y,z is 



/3Ni > 




8x 


lx 


9Z 




^ 9N 






3C 


’ 9C ’ 


3? 




9x 


3N . 
S 3n 




3x 

3n 


IX 

’ 3n ’ 


9z 

3n 




9N^. 

3y 


3N. 




3x 


3y 


9z 




9N . 
1 


J 




34 


’ 94 ’ 


34 




[3Z J 
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The Jacobian matrix [J] is defined as 



[J] 



3x 


lx 


3z 


35 


’ 3C ’ 


35 


3x 


lx 


3z 


3n 


’ 3n ’ 


3n 


3x 


lx 


3z 


3? 


’ 3? ’ 


3? 



( 49 ) 



From equation (47), the Jacobian matrix becomes 



[J] = 



9N. 3N. SN. 

? ^i ’ ? 9T^ -^i ’ • ^ ^i 



3N. 3N. 

E I X E 1- 

^ 3 n ""i ’ \ 3n 



3N . 



3N. 3N. 

? JT~ ^i ’ ? ^ 



3N 

i ’ ? ^ 



(50) 



The determinant of the Jacobian is used to transform the vol- 
ume of element in cartesian coordinates to local coordinates. 
For a volume of element [11], 



dxdydz = det [J] d^dnd? 



(51) 



Note that the determinant of the Jacobian is a variable for 
elements of curved geometry. Only in the case of straight- 
sided elements is the determinant of the Jacobian a constant. 

In the plane of the triangle, the area coordinates (L-]* 
12 * 13 ) number one more than the cartesian coordinates (x,y). 
Thus, is defined as a dependent variable. This establishes 
the origin of the Cn coordinate system at corner point 3, as 
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illustrated in figure 11. Recall that the Cn axes need 
not be orthogonal. As such 



? = 



n = 



L 

L 



1 

2 



Using equat i on ( 42 ) , 



L3 1 - C - n 

Applying the chain rule yields 

3N. 3N^. 3L.J 3N^. 31^ 3N^. 31^ 

W~ ^ W~ w~ 

3N. 3N. 3L, 3N. 3L. 3N . 3L. 

I I ^ ^ I o 

3 n 3L-| 3ri ^^2 



(52a) 

(52b) 



(52c) 



(53a) 



(53b) 



Using equations (52a), (52b), 
and (53b) gives 



3N . 
1 


3N, 


9N, 




9L, 


3 L 3 


3N . 


3N . 


3N . 


1 


1 


3ti 


- Tq - 


3 L 3 



and (52c) in equations (53a) 

(54a) 

(54b) 
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Figure 11. n? coordinates in a triangle 
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Now the Jacobian matrix can be evaluated using the shape func- 
tions N^. = N^. (L^ ,l 2 ,?) . Inserting equations (54a) and 

(54b) in equation (50) produces 







3N. 


3N. 




3N. 


3N. 




3N. 


3N. 






i 


34 


>4 


i ^*-1 




)yi 


i ^'-1 


54 ) 4 - 






3N. 


3N. 




3N. 


3N. 




3N. 


3N. 


[J(L^L2 


, 13 , 4 )] = 


i ^‘-2 




)x, 


i ^‘-2 




)^i 


’ ^^3T~ ~ 
i ^‘■2 








3N. 






3N. 






3N. 










9 






)y^- , 






i 




















(55) 


To 


summar i ze 


, suppose 


i t 


is required 


to 


transform the 



i ntegral 

I = J F' (L^ ,L 2 . 1 - 3>0 dxdydz (56) 

Vol 

to an integral entirely in terms of local coordinates. The 
determinant of equation (55) is then utilized to give 

' F'(L^, 12 , 13 , ?)det[J(L^, 12 , 13, ?)]dL^dL2dc 
■''00 (57) 

Equation (56) is now in a form suitable for numerical inte- 
gration. 

The development of coordinate transformation to this 
point has been under the general assumption that all the 
sides of the element are curved. In this work, the only 
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curved side is in the plane of the triangle. The sides of 
the element along the prism axis are straight. As such, a 
modified Jacobian matrix can be employed, thereby reducing 
the number of calculations to be performed. This modified 
Jacobian is the 2x2 matrix defined by 



[J 




9N . 
i:(- ■' 



i ^1 



2 (- 



3N . 



9N 



, 3L2 



3N . 
i ^'-1 



3N 

917^) Xi 

^^3 ^ i 



3N . 
i:(- ' 



3L 



1 



3N 
3N . 

3L7)^i 



(58) 



Along the prism axis or ^ direction, it can be shown that 

dz = ^ dc (59) 

where h = height of the element. The volume relationship 
is then given by 

dxdydz = | det[J*] dL^dL 2 ? (60) 

The integral of equation (56) then assumes the form of 



1 1 1-L^ 

■ •?/// F'(Li ,L2.L3,?)det[J*(L^ ,L2,L3,U]dL^ ,dL2,dC ( 61 ) 
-loo 

Equation (61) is the basis of numerical integration applied 
in this work. 

F. CONSTRUCTION OF ELEMENT MATRICES 

The system matrix operators can be constructed through 
the use of the global basis functions G, or through the 
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application of the element shape functions. Although the 
global basis functions were used to demonstrate the method 
of Galerkin, the construction of the system matrices in this 
work was achieved through element considerations. The solu- 
tion of the unknown variable ^ within the element domain 
was approximated by 

= Z N. i|;^(t) , i = l ,2,. . . ,n® (62) 

i = 1 

where n is the number of element nodal points, N^. are 
the element shape functions, and i|^^-(t) are the time-depen- 
dent nodal magnitudes of . The element contribution to 

the system matrix operators is defined by Galerkin's ortho- 
gonality condition expressed by 

/ N . R® dVol = 0 , j = l ,2, . . . ,n® (63) 

J 

Vol 

where R® is defined by replacing ^ with in equations 

(37) and (38), and the integration is over the element volume. 
Using equation (62), the element contribution is portrayed 
for the core by 

[ r N.N.dVol]{ii),®(t)}-[ / N.v2N.dVol]{(cl •/).}+[ / N .N.dVol]{(c2-4'®).} 

+ [ / N.N.Jlnd + N K^(t))dVol] {(c4-ij;®) .} 

+ [ f N.N.dVol] { (c5*i|^(f) ).dt'} = 0 (64) 

Vol 0 
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where i,j,k = 1,2,... ,15 , and cl, c2, etc., are constants 
at node i. The bracketed expressions represent square 
matrices of size 15x15, and the braced expressions represent 
column vectors of size 15x1. For the reflector, the last 
two terms of equation (64) are zero. Before any operation 
can be performed on equation (64), the nonlinear terms (last 
two terms) must be "linearized" and the term with the 
operator must be integrated by parts. 

1 5 

The nonlinear feedback term, 5,n[l + -1—2 ’ 

^ 0 k=i 

was linearized by using predicted values 

of the unknown function at time t. These predicted values 
come from the time integration scheme by Franke [7J. 
Basically, the integration scheme utilizes a predictor-cor- 
rector method which predicts values of the unknown function 
at the next time by using the derivatives of the function. 

p 

Adopting the predicted values enabled integration over 

space. The term in equation (64) involving the nonlinear 
feedback can be written as 

[ f NjN^.5,n(l + •••+ N.j ^^) )dVol ] {(c4*4^®)^.} 

Vol ° 



The last tern of equation (64) describes the delayed 



neutron contribution. In general, 
t 



-Xt 
e n 






+ e 



-'*n / 



Xt' , e,.,\,., 
e (t )dt 



(65) 



'n-1 



where n is number of time steps, t^ is the current time. 
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and IS the previous time. To approximate the integrals 

p 

in equation (65), the predicted values were employed 

in a simple trapezoidal rule. The trapezoidal rule was be- 
lieved to be sufficient since very small time steps were uti- 
lized in this work. For example, at time t 2 > 

(SUM^.) = (previous SUM^.)e~^^ + ^( ® ( t.| )+ 



where 

H = current time step = t 2 - t.j 
(SUM.) = e'^^2 ■/ ^ t|;.®(t')dt' 

1 I 

0 

1 ■ 1 

(previous SUM^) = e ^^1 J ® ( t ' ) dt ' 

0 

The previous sum is formed by accumulating the trapezoidal 
integration of each time step. In general then, for time t^ 

(SUM^.) = (previous SUM^.)e’^^+ ^ ( ^n_ 1 -) ^ 

The delayed neutron term in equation (64) can be expressed 
as 

[ / N.N.dVol] {(C5-SUM) .} 

Vol 
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I 







In order to bring equation (64) to final form, the 
operator was integrated by parts. Since the flux at the sur- 
face of the reactor is zero, integration by parts yields 

- 9N. 9N. 9N. 9N. 9N. 9N. 

NjV^N.dxdydz = - / (_ ^ H. ^ ^ ^ ^)dxdydz (67) 

Vol Vol 



In vector notation, equation (67) is expressed as 



J NjV^N^.dxdydz = 

Vol Vol 



/ 9N. 'i 

9x 



^ 9N. 9N 

r 1 1 . 

J ’ 9y" ’ 9z"^ ^ 



9N. 

1 , 



9N. 

i 

9y 

9N. 

J 

9z 



r dxdydz (68) 



Using the chain rule produces 



9 Ni >1 
9 )T 




9 x 
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’ 9 x 
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r 9 N. 9N X 
^917- 9 UJ^ 


9 N. 




9 L^ 

9 y 


9L2 
’ 9 y 






9 N, 9 N, 

^91:7 - -^4 


j 3 y 


> - 


, 0 




9 N. 








’ 9 z 




9 N. 


1 

1 92 i 
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, 0 




‘h) J 



(69) 



Letting the 3x3 matrix above be [B‘], the term becomes 



/ 

Vol 



N N^. dxdydz = 

9N. 9N. 9N, 

■/ "^^91:7 ■ 9Lr)’^9ur 
Vol ' ^ 



^ 9N. 9N. ' 



9N . 9N . T j 



9N. 

JL' 

■3 



9N. 

^9uJ ■ 

9N 

If 



f dxdydz 



(70) 
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where [B']^ is the transpose of [B‘] . By applying the 
chain rule in equation (47) and using equation (59), it can 
be shown that for this work 



[B-] = 



1 



1 



3 N, 9 N,t- ’ 3 N, 9 N,(. 

9T[^ ^15^ ^1 9L^ ^15^ 



1 



1 



9 N, 9 N,. ’ 9 N, 9 N,(. 



, 0 



, 0 



2 

’ h 



(71) 



0 , 0 

[B']^ can be derived from equation (71). 

Note from equation (64) that there are three basic ele- 
ment matrices which are defined as follows after applying 
equati on ( 60) : 



' f /// '’' N.N.det[J*]dL,dL2dc 



(72) 



-loo 



[GG 



. 1 1 1 -L, 9 N. 9 N. 9 N. 9 N . 9 N. , 

^317 ■ 3r7>><7r7 ■ 

- 1 0 0 



9L^ 9L3' ’ ^TIT^ ’ ' 3 ? 



9 N. 9 N. 

J 

TUT 3 L, ^ 



1 "^3 

9 N. 9 N. 

^ 9 L^ ■ 9L3 

9 N . 



) 



3? 



k det[J ]dL-|dL 2 dc 



(73) 






[GGGj^] = f /// NjN.in(l+ +...+N, 5 i(., 5 ))det[J JdL^dLjd; 



-loo 



( 74 ) 
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Element matrices ^G..] and [GG..] are independent of 

J * J 1 

time. However, [GGGj^.] is time dependent due to the utili- 
zation of the predicted values which changes with time. 

In terms of these three basic element matrices, equation (64) 
i s 



+ [GGGjj] {(c4-(i®)^.} + [G .] {{c5-SUM),l = 0 (75) 

The last two terms of equation (75) are zero for the reflector. 

G. CONSTRUCTION OF THE SYSTEM MATRICES 

The 15x15 coefficient element matrices were calculated 
according to equations (72), (73), and (74); and the results 
were collected element by element into the corresponding sys- 
tem coefficient matrices. The system coefficient matrix 

[BIGG] is developed from [G..] , [BIGGG] from [GG..] , 

J I J T 

and [BIGH] from [GGG..] . [BIGG] and [BIGGG] are inde- 

J * 

pendent of time and can be constructed once and for all from 
geometry considerations. [BIGH] is dependent on both geom- 
etry and time due to the time dependence of the predicted 
flux utilized in the feedback term. Thus, [BIGH] is recal- 
culated at each time increment. 

Non-zero contributions to a global nodal point I come 
only from adjacent elements sharing that same nodal point I. 
Thus, the system matrices are sparse and banded. The pro- 
cess of assembling contributions from element matrices 
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requires the identification of a local nodal point (i=l,2, 
...,15) with a global nodal point ( 1 = 1 ,2 , . . . NUMNP , where 
NUMNP is the total number of system nodes). This correspond- 
ence between element and global nodes is accomplished via 
the connectivity matrix. 

The formal treatment of the field equations in terms of 
the system coefficient matrices is described by the equation 



[BIGG] + [BIGGG] {(cl-ii;)^} + [BIGG] {(c2->|^)j} 

+ [BIGH] {(c4*ij;)j} + [BIGG] {(c5*SUM)j} = 0 (76) 



where the system matrices are NUMNP x NUMNP and the column 
vectors are of length NUMNP x 1. However, the direct appli- 
cation of equation (76) requires a large amount of computer 
storage. To take advantage of the sparsity of the system 
matrices, an optimum compacting scheme described by Ref. 8 
was employed. 

The concept behind OCS is simply to store only the non- 
zero terms of a coefficient matrix. OCS requires two integer 
arrays , say JB and NAME, and a vector of non-zero coeffi- 
cients of the square system matrix. For purposes of illus- 
tration, the square system matrix is called B, and the vector 

t h 

of non-zero coefficients of B is called BB. The i^ integer 
entry in the NUMNP x 1 JB vector is the number q^. . This 
number is defined by 



where 



i = l 

q. = 1 + 2 p. , i=l ,2,. . . , NUMNP (77) 

J = 1 

is the number of terms in the i^^ equation. In 
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r;f" 
















other words, is the number of nodes that the node 

"sees". JB is therefore a poin'fer vector of length NUMNP+1 
^ h 

whose i term locates the initial position in the BB vector 

of the contributing coefficients to the i^*^ equation. The 

NUMNP 

NAME vector of length Mxl , where M = Z p. , consists of 

i = l ^ 

NUMNP successive vector blocks of variable length p^. , 
i = 1 , 2 , . . . , NUMN P . The p^. integer numbers in the i^*^ block 

of NAME list the p^. contributors to the i equation. The 

Mxl BB vector contains the real non-zero coefficients of the 
NUMNPxNUMNP B matrix, arranged in the same contiguous block 

■f* h H 

arrangement as the NAME vector. The j term in the i 
block or BB(JB(I) + J-1) is B(I,K), where K = NAME( JB( I )+J-l ) . 
To illustrate, consider the grid shown in figure 12. The two 
array vectors and the coefficient vector matrix of non-zero 

<1 , 5, 10, 13, 18, 25, 30, 33, 38, 42> 

<1 ,2,4,5|2,3,1 ,6,5| — |9,8,6,5> 

<Bi 1 ,B^ 2 ’®1 4 ’^1 5 1 ^22*^23 ’^21 ’^26 ’^25 ^ ^ ®99 ’^98 ’^96 ’^95^ 

In this illustration, NUMNP = 9 and M = 41 [8]. 

In this work, a judicious method of numbering the system 
nodal points was adopted to further reduce computer storage 
requirements. Since the surface boundary nodes of the reac- 
tor represent zero neutron fluxes, the contributions of these 
nodes to interior or non-zero nodes can be discarded. Thus, 
only the interior nodal points need to be considered. These 
non-zero nodes were numbered first in the finite element mesh 



terms are; 
JB 

NAME = 
BB 
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Number - element number 



Figure 12. Sample grid used for illustrating OCS 
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used so that in the OCS, the number of non-zero nodes (NNZ) 
replaces NUMNP. 

The vectors of non-zero coefficients will be designated 
BIGG, BIGGG and BIGH since the square coefficient matrices 
described in equation (76) were not utilized. To illustrate 
the application of OCS in the system, the following sample 
program is given: 

DO 450 1 = 1 , NNZ 

JBB = JB(I) 

JE = JB( I+l )-l 
DY(I) = 0.0 
DO 500 J=JBB, JE 
LL = NAME(J) 

DY(I) = DY(I) + BIGG( J)*ij;(LL) + B I GGG ( J ) *cl ( LL ) * 
i^(LL) + BIGG( J)*c2(LL)*iJ;(LL) + B IGH ( J ) *c4 ( LL ) * 

4^(LL) + BIGG(J)*c5(LL)*SUM(LL) (78) 

500 CONTINUE 
450 CONTINUE 

where 

LL = nodal point "seen" by node I 

JE-JBB = total number of nodal points that node I "sees" 

DY(I) = summation of all contributions to node I and 

should sum to zero as stated by the field equa- 
tion 

The assumption of a homogeneous reflector was relaxed in the 
application of equation (78). Interface nodes were assigned 
properties of the core. Therefore, if LL is a reflector 
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node not on 
of equation 



the core reflector interface, the last two terms 
(78) are nonexistent. 
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IV. NUMERICAL INTEGRATION 



A. LINE AND AREA INTEGRATION 

The solution of the element matrix equations was achieved 
through numerical integration since an exact closed form solu 
tion cannot be established. The volume integration was accom 
plished by using a line integration in the ^-direction and an 
area integration in the plane of the triangle. The line inte 
gral is described by the Gaussian quadrature formula [12] 




1 

f(c)dc = 



-1 




Hkf(a^) 



(79) 



where n is the number of Gauss integration points 
= weighting coefficients 

f(aj^) = the function f(c) evaluated at Gauss point aj^ 

Table III lists aj^, and n [11]. 

The area integration was achieved by the equation 

m 

/ f(L^ ,L2.L3)dL^dL2 = E w^f ( L^ ^ L 2 "’ , L3^ (80) 

« 'i m=l 



where iii is the number of area integration points and w^ 
are the weights. The numerical integration points for the 
area integration are given in Table IV which was extracted 
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TABLE III 



Abscissae and Weight Coefficients 
of the Gaussian Quadrature Formula 



/ f(?)d; 

-1 



n 



E H^f(a,) 

k=l 





±o 




n 




2 




H 




0-57735 


02691 


89626 








1-00000 


00000 


00000 






* 


n 


■= 


3 








077459 


66692 


41483 








0-55555 


55555 


55556 


000000 


00000 


00000 


n 




4 


0-88888 


88888 


88889 


0 86113 


63115 


94053 








0-34785 


48451 


37454 


0 33998 


10435 


84856 








065214 


51548 


62546 


0 90617 


98459 


38664 


n 




5 


0-23692 


68850 


56189 


0-53846 


93101 


05683 








0-47862 


86704 


99366 


000000 


00000 


00000 


n 




6 


0-56888 


88888 


88889 


0-93246 


95142 


03152 








0-17132 


44923 


79170 


066120 


93864 


66265 








0-36076 


15730 


48139 


0 23861 


91860 


83197 


n 




7 


0-46791 


39345 


72691 


0 94910 


79123 


42759 








0-12948 


49661 


68870 


074153 


11855 


99394 








0-27970 


53914 


89277 


0-40584 


51513 


77397 








0-38183 


00505 


05119 


0-00000 


00000 


00000 


n 




8 


0-41795 


91836 


73469 - 


0 96028 


98564 


97536 








0 10122 


85362 


90376 


0-79666 


64774 


13627 








0-22238 


10344 


53374 


0-52553 


24099 


16329 








0-31370 


66458 


77887 


0-18343 


46424 


95650 


n 




9 


0-36268 


37833 


78362 


0 96816 


02395 


07626 








0 08127 


43883 


61574 


083603 


11073 


26636 








0-18064 


81606 


^4857 


0-61337 


14327 


00590 








0-26061 


06964 


02935 


032425 


34234 


03809 








0-31234 


70770 


40003 


000000 


00000 


00000 


n 


s 


10 


0-33023 


93550 


01260 


0-97390 


65285 


17172 








006667 


13443 


08688 


086506 


33666 


88985 








0-14945 


13491 


50581 


067940 


95682 


99024 








0-21908 


63625 


15982 


0 43339 


53941 


29247 








0 26926 


67193 


09996 


0 14887 


43389 


81631 








0-29552 


42247 


14753 
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TABLE IV 



Numerical Formulas for Triangles 



Order Fig. Error Points Triangular Weights 

Co-ordinates 2 





“ i. i. 0 
b o,Li 
r i.O.i 



i 

i 



Cubic 





Cubic 




Quintic 




R - m*) 



R = 



a 


i. i.i 




b 


i.i.O ^ 


1 


c 


o.J.i 


r 65 


d 


i.oj J 


1 


e 


1.0,6 'I 


1 


f 


0.1.0 


f A 


9 


0,0.1 J 


1 


a 


U.i 


0-225 


b 




'i 


c 


U..3..13 


t, 013239415 


d 


^1. a 


■1 


e 


.?2t Pz 


f 


Pt. l2. ^ 


Ij 012593918 


9 


^2* Pi' ^2 



with 

a, = 005971587 
fiy = 0 47014206 
ij * 0 79742699 
/?2 = 010I28651 
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from Ref, 11. The volume integration of the function 
f (L^ ,?) using numerical integration is therefore 




,L2,L2,?)dL-|dL2d? = 



n 



m 




m f ^1 ^ r ^ 1 

1 ’’-2 ’‘-3 ’ 



(81) 



where ~ manner of equation (81), the element 

matrices become 



where F and T are easily derived from equations (73) and 



evaluated at each integration point. They are not shown as 
such merely for the sake of convenience in writing the equa- 
tions, 

B. NUMBER OF INTEGRATION POINTS 

It is difficult to estimate the number of integration 
points required for good accuracy due to the complexity of 
the functions involved. The basic rule that the best number 
of integration points is found by trial and experience was 
adopted. For the [G..] element matrix, the numbers involved 




k=l m=l 



(74), respectively. Note that N^. and det[J ] are also 



J 1 
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I 



can be approximated. This was done by letting the determi- 
nant of the Jacobian equal to twice the area of the curved 
triangle (checking det[J*] at each integration point showed 
that this assumption was not too unreasonable). With det[J*] 
outside the integration process, [G..] can be solved in 

J * 

closed form by integrating out c from -1 to +1 and then apply- 
ing the closed form equation [ 12 ] 



where m.j , m^ , m^ are positive integer exponents and A is 
the area of the triangle. 

Five test points within the 15x15 element matrix were 
selected, as listed in Table V. The values obtained from 
the application of equation (85) are also given in Table V. 
Three sets of area integration points, each with a different 
number of c Gauss points, were used. These three sets of 
area integration points, given in Table IV, are: 

1. cubic order (4 points) 

2. cubic order (7 points) 

3. quintic order (7 points) 

Using each of the three integration points above with differ- 
ent ? Gauss points in equation (82) produced the results 
obtained in Table V. From these results, the quintic order 
area integration points were selected for the element matrix 

[G..] with three c Gauss points in the prism axes. This 
J ^ 

set of integration points was also used for [GGGj^-]. 




‘ L 2 L 3 ^ d(Area) 





( m.| +m2+m3+2 ) ! 



(85) 
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TABLE V. 



Selection of Integration Points for [G..] 



Area 



5 points 


points 


G(l,l) 


G(2,2) 


_ G(1,9) 


G(6,9) 


G(9,9) 


13 


cubic 
(4 pts) 


1415 


5278 


-2756 


5072 


8536 


5 


II 


1817 


5278 


-3170 


5072 


10240 


7 


II 


1817 


5278 


-3170 


5072 


10240 


3 


cubic 
(7 pts) 


3248 


8247 


-3114 


4960 


10243 


5 


II 


3249 


8247 


-3114 


4960 


10240 


7 


II 


3249 


8247 


-3114 


4960 


10240 


3 


quintic 
(7 pts) 


2464 


6600 


-3144 


5020 


10240 


5 


II 


2464 


6600 


-3144 


5020 


10240 


7 


II 


2464 


6600 


-3144 


5020 


10240 


Approximated values 


2500 


6700 


-3142 


5026 


10053 
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The numbers for the [G6..] element matrix could not 

J * 

be approximated. The best that could be done was to obtain 

an idea of the order of magnitude of this element matrix. 

Towards this end, a linear triangular element v/as assumed. 

Using linear approximation, the order of magnitude was found 
3 

to be about 10 . Using the three different area integration 
points mentioned above with varying ? Gauss points in equa 
tion (83) yielded the results given in Table VI. Further 
checks on the [GG..] element matrix showed that the cubic 

<J * 

order with four area integration points yielded the desired 

order of magnitude. Thus, the fourth order cubic with five 

C Gauss points was employed for the [GG..] element matrix 

J ^ 

A note should be mentioned here in regards to the vast dif- 
ference in results obtained for [GG^^.] using different 
area integration points. Most likely, it was due to the 

[B'] and [B’]'*' matrices which required the inversion of 
3x 



3x 

3U 



3L, 



, etc. Or perhaps it was caused by the nature of 



1 

the hybrid element used in this work. In any case, further 
investigation is warranted in this area. 
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;.J. . 




TABLE VI 



Selection of Integration Points for [26^^.] 



Area 



C, 1 II uo 

3 


1 II uo 

cubic 
(4 pts) 


\ 1 j -/ y 

18.4 


uu V ^ 5 o y 

170.1 


UU V 1 •I’ 5 j y 

-165 


UU V u 5 u y 

165.9 


Cl 13 \ \ ^ J 

106.1 


5 


II 


22.0 


177.6 


-186 


195.9 


106.1 


7 


II 


22.0 


177.6 


-186 


195.9 


106.1 


3 


cubic 
(7 pts) 


-1.03x10^ 


2.11x10^ 


1.05x10^ 


8.42x10^ 


2056 


5 


M 


-1.03x10^ 


2.11x10^ 


1.05x10^ 


8.42x10^ 


2056 


7 


11 


-1.03x10^ 


2.11x10^ 


1.05x10^ 


8.41x10^ 


2056 


3 


quintic 
(7 pts) 


-52.6 


1.27x10^ 


-6.64x10^ 


6.84x10^ 


1.51x10^ 


5 


11 


-52.6 


1.27x10^ 


-6.64x10^ 


6.84x10^ 


1.51x10^ 


7 


II 


-52.6 


1.27x10^ 


-6.64x10^ 


6.84x10^ 


1.51x10^ 
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V. TEST PROBLEMS AND RESULTS 



The reactor was subjected to uniform and local perturba- 
tions in the form of a ramp input described by 

+ at (86) 

★ 

where a is the change in Z^ per unit time and Z^ is 

★ 

the critical fission cross-section. Z^ must first be ob- 
tained before the perturbations can be applied. This was ac 

complished by trial and error until a stationary solution 

* 

was reached. For mesh I, Z^ was found to be 0.0057360 

★ 

per cm. No attempt was made to find Z^ for mesh II due 
to time limitations. As such, the test problems outlined be 
low were applied to mesh I. Future work is planned to apply 
the test problems on mesh II. 

The following perturbati ons were applied: 

a) Uniform perturbation of 10 dollar of reactivity per 
second : 

Z^(r,t) = Z^ + at , in the core 

where a = 0 . 005893/cm-sec 

b) Local perturbation at the core center of 100 dollar 
of reactivity per second: 

Z^(r,t) = Z^ +at6(r^) , in the core 
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where 



is (0,0,0) and a = 0.015407/ cm -sec 



c) Local, off-center perturbation: 

2f(r,t) = + a^.t6(r.|) , i = 1,2,3, in the core 

where 

r^ = (0,60,40) 



a.| = 0 . 01 5407/ cm- sec = 


100 


dollar 


per 


second 


a 2 = 0 . 0081 23/ cm- sec = 


50 


dollar 


per 


second 


- 0 . 005894/cm-sec = 


10 


dollar 


per 


second 



Three test points, (0,0,0), (60,0,0), and (-60,0,80), 
were selected to trace the neutron time history. For cases 
a) and b), the neutron flux was plotted at each test point 
during transience. This is shown in figures 13 and 14. 

Case c) involved three ramp inputs and were conducted for 
both the linear and nonlinear reactor equations. The linear 
and nonlinear responses were compared at each test point for 
each ramp input and are illustrated in figures 15 through 23. 
The radial and axial flux distributions at time t = 0.0123 
second were plotted for the steady state and the 100 dollar 
perturbation case. These are embodied in figures 24 and 25. 
Finally, a neutron flux early time history between mesh I 
and mesh II was plotted to check the effect of using a finer 
element mesh. The result is portrayed in figure 26. 

Figures 13 and 14 revealed a clear space dependence of 
the neutron flux during transience as expected. Figures 15 
through 23 demonstrated the effects of temperature feedback 
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neutron flux 




time (sec) 



A - neutron flux at test point (0,0,0) 

B - neutron flux at test point (60,0,0) 

C - neutron flux at test point (-60,0,80) 



Figure 13. Neutron flux transient history at three 
test points with a uniform perturbation 
of 10 dollar of reactivity per second. 
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neutron flux 




A - neutron flux at test point (0,0,0) 

B - neutron flux at test point (60,0,0) 

C - neutron flux at test point (-60,0,80) 



Figure 14. Neutron flux transient history at various 

test points for a local central perturbation 
of 100 dollar/sec of reactivity. 
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neutron flux 




A - nonlinear neutron response 
B - linear neutron response 



Figure 15. Linear and nonlinear fluxes at (0,0,0) due 

to a local perturbation of 100 dollar/sec of 
reactivity at (0,60,40). 
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neutron flux 




time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 16. Linear and nonlinear fluxes at (60,0,0) due 
to a local perturbation of 100 dollar/sec at 
(0,60,40) . 
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neutron flux 




lo“^ lo"'^ lo"^ lo' 

time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 17. Linear and nonlinear fluxes at (-60,0,80) 

due to a local perturbation of 100 dollar/sec 
at (0,60,40). 
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10 



9 



—I 1 1 1. 

-5 -4 -3 

10 10 10 10 

time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 18. Linear and nonlinear fluxes at (0,0,0) due 
to a 50 dollar/sec local perturbation at 
(0,60,40) . 
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neutron flux 




10~^ lo"'^ 10~^ 10 ^ 



time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 19. Linear and nonlinear fluxes at (60,0»0) 

due to a 50 dollar/sec local perturbation 
at (0,60,40). 
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neutron flux 




time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 20. Linear and nonlinear fluxes at (-60,0,80) 
due to a 50 dollar/sec local perturbation 
at (0,60,40). 
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neutron flux 




time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 21. Linear and nonlinear fluxes at (0,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40). 
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neutron flux 




A - nonlinear neutron response 
B - linear neutron response 



Figure 22. Linear and nonlinear fluxes at (60,0,0) due 
to a 10 dollar/sec local perturbation at 
(0,60,40) . 
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neutron flux 




-5 -4 -3 

10 10 10 10 

time (sec) 

A - nonlinear neutron response 
B - linear neutron response 



Figure 23. Linear and nonlinear fluxes at (-60,0,80) 
due to a 10 dollar/sec local perturbation 
at (0,60,40). 
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neutron flux 




I — I 1 (— ► 



60 30 0 -30 -60 

x=0, z=0, y=y 

# - point of local perturbation 
A - neutron flux at steady-state 
B - neutron flux under local perturbation 

Figure 24. Radial flux distribution for the steady state 
and 100 dollar/sec local perturbation. 
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neutron flux 




< 1— H 1 1 1 1 ( 

95 80 40 0 -40 -80 -95 

X=0, y=60, Z=Z 

0 - point of local perturbation 

A - axial flux distribution during steady state 
B - axial flux distribution during perturbation 

Figure 25. Axial flux distribution for the steady state and 
100 dollar/ sec local perturbation. 
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X 




time (sec) 



A - mesh I 
B - mesh II 



Figure 26. Early time history of the neutron flux at core 
center using mesh I and mesh II. 
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>1 



and delayed neutrons on the flux. To obtain an idea of the 

difference in magnitude between the linear and nonlinear 

flux for the 100 dollar local perturbation, at time t = 10"^ 

second, the linear case predicted a neutron flux at the core 

12 2 

center of 2.024 x 10 neutrons per cm per sec. The non- 

10 2 

linear reactor predicted 1.85 x 10 neutrons per cm per 
sec. Although the numbers are not sufficiently accurate due 
to the crudeness of the finite element mesh, the significant 
difference in the order of magnitude between the linear and 
nonlinear neutron flux leads to the belief that the three- 
dimensional finite element model utilized is predicting the 
correct trend of neutron flux behavior. 

The radial flux distribution for the steady state shown 
in figure 24 portrays the expected flux distribution. The 
radial flux distributions for the local perturbation case 
show a skew distribution at the point of perturbati on . The 
axial flux distribution at steady state is very much symme- 
tric about the center. Again, the expected skew at the 

point of perturbation is there, as shown in figure 25. 

★ 

Figure 26 gives an indication that the for finer 

meshes is higher. However, curve B of figure 26 should be 
extended to longer times to verify this hypothesis. Curve B, 
as plotted in figure 26, used four hours of computer time 
employing the H-compiler of the IBM 360/67. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



The finite element mesh employed here is crude, and, 
thus, results that were obtained should be considered as posi- 
tive indicators rather than numerically conclusive facts. The 
trend of neutron flux behavior is the major thrust of this 
work. From the results obtained, it is concluded that the 
expected patterns of neutron behavior as predicted by the 
three-dimensional finite element model used do occur. These 
patterns were best demonstrated by the differences between 
the linear and nonlinear flux responses and by the spatial 
flux distribution at steady state and during local perturba- 
tion. The three-dimensional quadratic finite element model 
utilized should produce better results by resorting to a 
finer mesh. The draw-back to a finer mesh, of course, is 
the significant increase in computer time and storage re- 
quirements, Once accurate results are obtained through 
finer element meshes, comparisons between three- and two- 
dimensional models can be attempted, 

A mesh generator was not developed for this problem. It 
is recommended that this be done to ease the transition from 
one mesh to another and to minimize human error. In addi- 
tion, a similar calculation using a three-dimensional linear 
element should be performed to corroborate the results ob- 
tained here. It should be particularly noted that this type 
of problem is highly sensitive to the fission cross-section. 
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★ 

In the search for , it was necessary to adjust to 

the sixth decimal place. There is no exact method of deriv- 
★ 

ing due to the highly nonlinear aspect of the problem; 

therefore, trial and error must be used. Finally, the 
Gaussian quadrature used to determine the element matrices 
should be investigated. The cause of the differences in 
values obtained by using different number of integration 
points should be established. Was it due to the integration 
points, the coordinate transformation, or the element shape 
functions? This is an important question upon which future 
numerical results could be based. 
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The coordinates of the 128-element mesh are as follows 
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